# Setup -------------------------------------------------------------------
rm(list = ls())
source("figure_design.R")

# Data --------------------------------------------------------------------

vdem <- readRDS("vdem_v14_subset.rds")
covariates <- readRDS("covariates_subset.rds")
gallup <- readRDS("gallup_subset.rds")

# Code V-Dem data ---------------------------------------------------------

vdem_subset <- vdem %>%
  filter(year>2008) %>%
  left_join(covariates, by = c("country_name","year")) %>%
  mutate(v2x_regime = case_when(
    v2x_regime == 0 ~ "Closed Autocracy" ,
    v2x_regime == 1 ~ "Electoral Autocracy",
    v2x_regime == 2 ~ "Electoral Democracy",
    v2x_regime == 3 ~ "Liberal Democracy",
    TRUE ~ NA_character_),
    country_name = case_when(country_name == "Czechia" ~ "Czech Republic",
                             country_name == "United Kingdom" ~ "United Kingdom",
                             country_name == "United States of America" ~ "United States",
                             country_name == "Burma/Myanmar" ~ "Myanmar",
                             country_name == "Sao Tome and Principe" ~ "Sao Tome & Principe",
                             country_name == "Türkiye" ~ "Turkey",
                             TRUE ~country_name)) %>%
  rename(country = country_name,
         regime = v2x_regime,
         democracy_score = v2x_polyarchy) %>%
  select(country,year,regime,democracy_score,consumer_price_index,unemployment,gdp,gdp_growth)

vdem_subset %>% select(year, country, regime) %>% 
  filter(regime!="Closed Autocracy") %>%
  filter(year>2008) %>%
  mutate(state_2011 = case_when(regime=="Electoral Autocracy"&year==2011~1,
                                (regime=="Electoral Democracy"|regime=="Liberal Democracy")&year==2011~0,
                                TRUE ~ NaN),
         state_2020 = case_when(regime=="Electoral Autocracy"&year==2020~1,
                                (regime=="Electoral Democracy"|regime=="Liberal Democracy")&year==2020~0,
                                TRUE ~ NaN)) %>%
  filter(!(is.na(state_2020)&year==2020)) %>%
  group_by(country) %>%
  summarise(status_2011 = sum(state_2011,na.rm=T),
            status_2020 = sum(state_2020,na.rm=T)) %>%
  mutate(long_regime = case_when(status_2011==1&status_2020==1~"Stable Autocracy",
                                 status_2011==1&status_2020==0~"Democratization",
                                 status_2011==0&status_2020==0~"Stable Democracy",
                                 status_2011==0&status_2020==1~"Autocratization"))-> 
  vdem_regime_universe

vdem_regime_universe %>% filter(long_regime=="Stable Autocracy") -> autocracy
vdem_regime_universe %>% filter(long_regime=="Stable Democracy") -> democracy
vdem_regime_universe %>% filter(long_regime=="Autocratization") -> autocratize
vdem_regime_universe %>% filter(long_regime=="Democratization") -> democratize

list_study_universe <- levels(as.factor(vdem_subset$country))
list_autocracy <- levels(as.factor(autocracy$country))
list_democracy <- levels(as.factor(democracy$country))
list_democratize <- levels(as.factor(democratize$country))
list_autocratize <- levels(as.factor(autocratize$country))

vdem_subset %>%
  left_join(vdem_regime_universe, by = "country") ->
  vdem_subset_complete

# Recode survey data ------------------------------------------------------

remove(vdem_regime_universe,vdem_subset,autocracy,autocratize,democracy,democratize,vdem,covariates)

gallup %>%  
  rename(country = WP5,
         year = YEAR_WAVE,
         trust_elections = WP144,
         leader_approval = WP13125) %>%
  mutate(country = remove_labels(to_factor(country)),
         year = remove_labels(year)) %>%
  mutate(trust_elections=case_when(trust_elections==1~1,
                                   trust_elections==2~0,
                                   TRUE ~ NA_integer_),
         leader_approval=case_when(leader_approval==1~1,
                                   leader_approval==2~0,
                                   TRUE ~ NA_integer_),
         leader_approval_na=case_when(leader_approval==1~1,
                                  TRUE ~ 0)) %>%
  filter(!is.na(leader_approval)) %>%
  mutate(leader_approval_bin = case_when(leader_approval == 1 ~ "Yes",
                                         leader_approval == 0 ~ "No"),
         leader_approval_bin_na = case_when(leader_approval == 1 ~ "Yes",
                                            leader_approval == 0 ~ "No"),
         leader_approval_bin = factor(leader_approval_bin, levels = c("Yes","No")),
         leader_approval_bin_na = factor(leader_approval_bin_na, levels = c("Yes","No"))) %>%
  left_join(vdem_subset_complete,by=c("country"="country","year"="year")) %>%
  mutate_at(vars(c(consumer_price_index,unemployment,gdp,gdp_growth)), ~ (. - mean(.,na.rm=T)) / (2 * sd(.,na.rm=T))) %>%
  filter(country %in% list_study_universe) ->
  regression_data

# Figure 3a ---------------------------------------------------------------

regression_data  %>%
  filter(year > 2011) %>%
  drop_na(long_regime,leader_approval,country,year) %>%
  group_by(country,year) %>%
  summarise(country = first(country),
            year = first(year)) %>%
  summarise(n())

regression_data  %>%
  filter(year > 2011) %>%
  drop_na(long_regime,leader_approval) %>%
  as_survey(weights = wgt) %>%
  group_by(year,long_regime,leader_approval)  %>%  
  summarise(mean=survey_mean(trust_elections,na.rm=T,vartype=c("ci")))  %>% 
  mutate(leader_approval = case_when(leader_approval==1~"Support Leader",
                                     leader_approval==0~"Oppose Leader"),
         leader_approval = factor(leader_approval, levels = c("Support Leader","Oppose Leader") )) -> 
  mean_data

plot_by_supporters <- ggplot(mean_data, aes(x=year, y=mean,color=leader_approval)) +
  geom_line(aes(linetype = leader_approval),size=0.8) + 
  geom_point() + xlab("Year") +
  ylab("Confidence in Elections") +  
  theme(legend.title=element_blank()) +
  geom_errorbar(aes(ymin=mean_low, ymax=mean_upp), width=.05)  + 
  scale_x_continuous(breaks = c(2011,2013,2015,2017,2019,2021)) + 
  facet_grid(~long_regime) + 
  theme(legend.position="top")  + 
  scale_color_manual(values = plot_colors)

plot_by_supporters

ggsave(plot_by_supporters, filename = "fig_3a.pdf", device = cairo_pdf, 
       width = 8, height = 3, units = "in")

# Figure 3a ---------------------------------------------------------------

mean_data %>% select(year, long_regime, leader_approval, mean) %>%
  pivot_wider(names_from = leader_approval, values_from = mean) %>%
  mutate(differential = `Support Leader` - `Oppose Leader`) ->
  mean_data_differential

differential <- ggplot(mean_data_differential, aes(x=year, y=differential)) +
  geom_line(size=0.8) + 
  geom_point() + xlab("Year") +
  ylab("Confidence Differential") +  
  theme(legend.title=element_blank()) +
  scale_x_continuous(breaks = c(2011,2013,2015,2017,2019,2021)) + 
  facet_grid(~long_regime) + 
  theme(legend.position="top")  + 
  scale_color_manual(values = plot_colors)

ggsave(differential, filename = "fig_3b.pdf", device = cairo_pdf, 
       width = 8, height = 2.5, units = "in")

# Regression models -------------------------------------------------------

regression_data %>%
  mutate(year = as.factor(year),
         country = as.factor(country)) %>%
  drop_na(trust_elections,leader_approval_bin,democracy_score,consumer_price_index,unemployment,gdp,gdp_growth,wgt)->
  regression_data_2

regression_data_2 %>%
  group_by(country,year) %>%
  summarise(country = first(country),
            year = first(year)) %>%
  summarise(n()) ->
  n_countries

summary(n_countries$`n()`)

regression_data_2 %>%
  group_by(country,year) %>%
  summarise(n()) ->
  n_sample

summary(n_sample$`n()`)

regression_data_2 %>%
  summarize(num_countries = n_distinct(country)) %>%
  pull(num_countries) ->
  n_countries

m1 <- glm(trust_elections ~ leader_approval_bin +
            consumer_price_index + unemployment + gdp + gdp_growth + 
            country + year, 
          weights = wgt,
          family = binomial(link = "probit"), data = regression_data_2)

# Compute clustered standard errors for model m1
m1_clust_se <- coeftest(m1, vcov = vcovCL(m1, cluster = regression_data_2$country))
# Replace standard errors in model m1 with clustered standard errors
m1$se <- m1_clust_se[, "Std. Error"]

m2 <- glm(trust_elections ~ leader_approval_bin + democracy_score +
            democracy_score/leader_approval_bin +
            consumer_price_index + unemployment + gdp + gdp_growth +
            country + year, 
          weights = wgt, family = binomial(link = "probit"), data = regression_data_2)

# Compute clustered standard errors for model m2
m2_clust_se <- coeftest(m2, vcov = vcovCL(m2, cluster = regression_data_2$country))
# Replace standard errors in model m2 with clustered standard errors
m2$se <- m2_clust_se[, "Std. Error"]

# Table 3 -----------------------------------------------------------------

texreg(list(m1,m2),
      file = "tab_3.tex",
       include.ci = F,
       digits = 3,
       custom.coef.map = list(
         "leader_approval_binNo" = "Oppose leader (vs. support leader)",
         "democracy_score" = "Polyarchy",
         "leader_approval_binNo:democracy_score" = "Oppose leader $\\times$ polyarchy",
         "consumer_price_index" = "Consumer price index",
         "unemployment" = "Unemployment",
         "gdp" = "Log(GDP)",
         "gdp_growth" = "GDP growth"
       ),
       groups = list("Country-level" = 4:7),
       booktabs = T,
       use.packages=F,	
       custom.gof.rows = list("Country Fixed Effects" = c("$\\checkmark$","$\\checkmark$"),
                              "Year Fixed Effects" = c("$\\checkmark$","$\\checkmark$"),
                              "$N$ Countries" = c(n_countries,n_countries)),
       include.rmse = FALSE,
       caption = "Two-way fixed effects probit regression: predicting confidence in elections with country’s electoral democracy score and individual leader approval. Standard errors in parenthesis. Data source: \\textit{Gallup World Poll}.",
       label = "tab3")

# Figure 4a ---------------------------------------------------------------

m2_me <- ggeffects::ggpredict(m2,c("democracy_score[0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]","leader_approval_bin")) %>%
  as.data.frame() %>%
  mutate(group = case_when(group == "Yes" ~ "Support Leader",
                           group == "No" ~ "Oppose Leader"),
         group = factor(group, levels = c("Support Leader", "Oppose Leader")))

plot_me <- ggplot(m2_me,aes(y=predicted,x=x,shape=group, fill=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) + 
  geom_line(aes(linetype = group, color = group)) +
  labs(x="Polyarchy", y="Pr(Confidence in Elections)") +
  theme(legend.position="top") +
  scale_color_manual(values = plot_colors) +
  scale_fill_manual(values = plot_colors)

plot_me

ggsave(plot_me, filename = "fig_4a.pdf", device = cairo_pdf, 
       width =  5, height = 4, units = "in")

# Calculate the percentage change in predicted values for "Oppose Leader"
# Identify the maximum predicted value at Polyarchy score = 1
opp_max <- m2_me %>% filter(x == 1 & group == "Oppose Leader") %>% pull(predicted)
# Identify the minimum predicted value at Polyarchy score = 0.4
opp_min <- m2_me %>% filter(x == 0.4 & group == "Oppose Leader") %>% pull(predicted)

# Percentage point change
round((opp_min - opp_max) * 100,1)
# Calculate the percentage change from opp_max to opp_min for "Oppose Leader"
round((((opp_min - opp_max) / opp_max) * 100), 1)

# Repeat the percentage change calculation for "Support Leader"

# Identify the maximum predicted value at Polyarchy score = 1
opp_max <- m2_me %>% filter(x == 1 & group == "Support Leader") %>% pull(predicted)
# Identify the minimum predicted value at Polyarchy score = 0.4
opp_min <- m2_me %>% filter(x == 0.4 & group == "Support Leader") %>% pull(predicted)

# Percentage point change
round((opp_min - opp_max) * 100, 1)
# Calculate the percentage change from opp_max to opp_min for "Support Leader"
round((((opp_min - opp_max) / opp_max) * 100), 1)

# Figure 4b ---------------------------------------------------------------

# Define the conditional effects function for probit models
conditional.effects_probit <- function(model, x, m, quantiles = 10) {
  interact <- paste0(x, ':', m)
  beta.hat <- coef(model)
  covs <- vcov(model)
  
  # Compute the range of values for the moderator variable
  z0 <- quantile(model$model[, m], seq(0, 1, 1 / quantiles))
  
  # Compute the conditional effects
  dy.dx <- beta.hat[x] + beta.hat[interact] * z0
  se.dy.dx <- sqrt(covs[x, x] + z0^2 * covs[interact, interact] + 2 * z0 * covs[x, interact])
  
  # Compute the predicted probabilities at each value of the moderator variable
  probs <- pnorm(dy.dx / sqrt(1 + se.dy.dx^2))
  
  # Compute confidence intervals using the standard normal distribution
  p_se <- dnorm(dy.dx / sqrt(1 + se.dy.dx^2)) * se.dy.dx
  upr <- probs + 1.96 * p_se
  lwr <- probs - 1.96 * p_se
  
  data.frame(m = z0, b = probs, lwr, upr)
}

con_effects <- conditional.effects_probit(m2, 
                                          x = "leader_approval_binNo", 
                                          m = "democracy_score", 
                                          quantiles = 500) %>%
  as.tibble

scatter_plot <- ggplot(con_effects, aes(x = m, y = b)) +
  geom_point(color = "white") +
  geom_smooth(aes(ymin = lwr, ymax = upr), stat = "identity", color = "black", linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  labs(x = "Polyarchy", y = "Effect of Oppose Leader")

# Add marginal histogram of x-axis variable
mod_plot_with_hist <- ggMarginal(scatter_plot,
                                 margins = "x",
                                 type = "density",
                                 color = "white",
                                 fill = "darkgrey",
                                 size = 3) 

mod_plot_with_hist
dev.copy(pdf, "fig_4b.pdf",
         width =  5, height = 4)
dev.off()

# APPENDIX ----------------------------------------------------------------

# Table B.3 ---------------------------------------------------------------

regression_data_2 %>% select(trust_elections,leader_approval,
                  consumer_price_index,unemployment,gdp,gdp_growth) %>% 
  as.data.frame() ->
  descriptives

stargazer(descriptives,type = "latex",
          digits.extra = 3,
          out = "tab_b3.tex",
          covariate.labels = c("Trust elections","Leader approval","Consumer price index","Unemployment","Log(GDP)","GDP growth"))


# Table B.4 ---------------------------------------------------------------

regression_data_2 %>%
  group_by(country,long_regime) %>%
  summarise(n = n()) ->
  list_countries

write.csv(list_countries, file = "tab_b4.csv", row.names = F)

# Figure B.1 --------------------------------------------------------------

regression_data_2 %>% select(country,year,trust_elections,democracy_score,
                           leader_approval_bin,consumer_price_index,unemployment,
                           gdp,gdp_growth) %>% 
  mutate(leader_approval_bin = case_when(leader_approval_bin == "Yes" ~ 0,
                                         leader_approval_bin == "No" ~ 1,
                                         TRUE ~ NA_real_)) %>% 
  drop_na(trust_elections,democracy_score,leader_approval_bin,consumer_price_index,unemployment,
          gdp,gdp_growth) %>% 
  as.data.frame() -> 
  binning_data

binning <- interflex(Y = "trust_elections",
                     D = "leader_approval_bin", 
                     X = "democracy_score", 
                     data = binning_data,
                     estimator = "binning",
                     na.rm = TRUE,
                     Xlabel = "Polyarchy",
                     Ylabel = "Confidence in elections",
                     Dlabel = "Oppose leader",
                     FE = c("country", "year"), 
                     Z = c("consumer_price_index","unemployment","gdp",
                           "gdp_growth"),
                     theme.bw = T)

ggsave(plot(binning), filename = "fig_b1.pdf", device = cairo_pdf, 
       width = 5, height = 7, units = "in")

# Table B.5 ---------------------------------------------------------------

regression_data %>%
  mutate(year = as.factor(year),
         country = as.factor(country)) %>%
  drop_na(trust_elections,leader_approval_bin_na,
          democracy_score,consumer_price_index,unemployment,gdp,gdp_growth,wgt)->
  regression_data_3

m1 <- glm(trust_elections ~ leader_approval_bin_na +
            consumer_price_index + unemployment + gdp + gdp_growth + 
            country + year, 
          weights = wgt,
          family = binomial(link = "probit"), data = regression_data_3)

# Compute clustered standard errors for model m1
m1_clust_se <- lmtest::coeftest(m1, vcov = vcovCL(m1, cluster = regression_data_3$country))
# Replace standard errors in model m1 with clustered standard errors
m1$se <- m1_clust_se[, "Std. Error"]

m2 <- glm(trust_elections ~ leader_approval_bin_na + democracy_score +
            democracy_score/leader_approval_bin_na +
            consumer_price_index + unemployment + gdp + gdp_growth +
            country + year, 
          weights = wgt, family = binomial(link = "probit"), data = regression_data_3)

# Compute clustered standard errors for model m2
m2_clust_se <- lmtest::coeftest(m2, vcov = vcovCL(m2, cluster = regression_data_3$country))
# Replace standard errors in model m2 with clustered standard errors
m2$se <- m2_clust_se[, "Std. Error"]

texreg(list(m1,m2),
          file = "tab_b5.tex",
          include.ci = F,
          digits = 3,
          custom.coef.map = list(
            "leader_approval_bin_naNo" = "Oppose leader (vs. support leader)",
            "democracy_score" = "Polyarchy",
            "leader_approval_bin_naNo:democracy_score" = "Oppose leader $\\times$ polyarchy",
            "consumer_price_index" = "Consumer price index",
            "unemployment" = "Unemployment",
            "gdp" = "Log(GDP)",
            "gdp_growth" = "GDP growth"
          ),
          groups = list("Country-level" = 4:7),
          booktabs = T,
          use.packages=F,	
          custom.gof.rows = list("Country Fixed Effects" = c("$\\checkmark$","$\\checkmark$"),
                                 "Year Fixed Effects" = c("$\\checkmark$","$\\checkmark$"),
                                 "$N$ Countries" = c(n_countries,n_countries)),
          include.rmse = FALSE,
          caption = "Two-way fixed effects probit regression: predicting confidence in elections with country’s electoral democracy score and individual leader approval. Standard errors in parenthesis. Data source: \\textit{Gallup World Poll}.",
          label = "tab3")

# Figure B.2 ---------------------------------------------------------------

m2_me <- ggeffects::ggpredict(m2,c("democracy_score[0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]","leader_approval_bin_na")) %>%
  as.data.frame() %>%
  mutate(group = case_when(group == "Yes" ~ "Support Leader",
                           group == "No" ~ "Oppose Leader"),
         group = factor(group, levels = c("Support Leader", "Oppose Leader")))

plot_me <- ggplot(m2_me,aes(y=predicted,x=x,shape=group, fill=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) + 
  geom_line(aes(linetype = group, color = group)) +
  labs(x="Polyarchy", y="Pr(Confidence in Elections)") +
  theme(legend.position="top") +
  scale_color_manual(values = plot_colors) +
  scale_fill_manual(values = plot_colors)

plot_me

ggsave(plot_me, filename = "fig_b2a.pdf", device = cairo_pdf, 
       width =  5, height = 4, units = "in")

# Figure 4b ---------------------------------------------------------------

# Define the conditional effects function for probit models
conditional.effects_probit <- function(model, x, m, quantiles = 10) {
  interact <- paste0(x, ':', m)
  beta.hat <- coef(model)
  covs <- vcov(model)
  
  # Compute the range of values for the moderator variable
  z0 <- quantile(model$model[, m], seq(0, 1, 1 / quantiles))
  
  # Compute the conditional effects
  dy.dx <- beta.hat[x] + beta.hat[interact] * z0
  se.dy.dx <- sqrt(covs[x, x] + z0^2 * covs[interact, interact] + 2 * z0 * covs[x, interact])
  
  # Compute the predicted probabilities at each value of the moderator variable
  probs <- pnorm(dy.dx / sqrt(1 + se.dy.dx^2))
  
  # Compute confidence intervals using the standard normal distribution
  p_se <- dnorm(dy.dx / sqrt(1 + se.dy.dx^2)) * se.dy.dx
  upr <- probs + 1.96 * p_se
  lwr <- probs - 1.96 * p_se
  
  data.frame(m = z0, b = probs, lwr, upr)
}

con_effects <- conditional.effects_probit(m2, 
                                          x = "leader_approval_bin_naNo", 
                                          m = "democracy_score", 
                                          quantiles = 500) %>%
  as.tibble

scatter_plot <- ggplot(con_effects, aes(x = m, y = b)) +
  geom_point(color = "white") +
  geom_smooth(aes(ymin = lwr, ymax = upr), stat = "identity", color = "black", linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  labs(x = "Polyarchy", y = "Effect of Oppose Leader")

# Add marginal histogram of x-axis variable
mod_plot_with_hist <- ggMarginal(scatter_plot,
                                 margins = "x",
                                 type = "density",
                                 color = "white",
                                 fill = "darkgrey",
                                 size = 3) 

mod_plot_with_hist
dev.copy(pdf, "fig_b2b.pdf",
         width =  5, height = 4)
dev.off()